# Import PAM data
um <- read_csv("data/20201014_UM_Acer/PAM/UM_pamdata.csv")
crf <- read_csv("data/20201013_CRF_Acer/PAM/CRF_pamdata.csv")
rrt <- read_csv("data/20201007_RRT_Acer/PAM/RRT_Acer_pamdata.csv")
mote <- read_csv("data/20201011_Mote_Acer/PAM/Mote_pamdata.csv")
fwc <- read_csv("data/20201010_FWC_Acer/PAM/FWC_pamdata.csv")
df <- list(um = um, crf = crf, rrt = rrt, mote = mote, fwc = fwc) %>%
bind_rows(.id = "nursery")
# Tidy PAM data - initial filtering
dff <- df %>%
# Make missing values zero (high temp / no signal)
replace_na(replace = list(fvfm = 0)) %>%
# Omit data for one coral that had missing data at 32°C
filter(!(max_temp == 32 & fvfm == 0)) %>%
# Change abnormally high values at high temps to zero
mutate(fvfm = if_else(max_temp > 35 & fvfm >= 0.4 & f <= 250, 0, fvfm)) %>%
# Discard all values > 0.750
filter(fvfm < 0.750)
# Plot all data points
ggplot(dff, aes(x = max_temp, y = fvfm)) +
geom_point(alpha = 0.2)
# Fit 3-parameter LL model to each coral individually (faster than all together) for initial outlier detection
# Define function to fit 3-parameter LL model to data and return NULL if fitting error
ll3 <- function(data) {
drm(fvfm ~ max_temp, data = data,
fct = LL.3(names = c("hill", "max", "ED50")),
upperl = c(100, 0.67, NA),
lowerl = c(NA, 0.66, NA))}
tryll3 <- possibly(ll3, otherwise = NULL)
# Fit model to each coral, get parameters, fitted values, and residuals
initmods <- dff %>%
nest(data = c(f, fm, fvfm, tank_id, tank_set, max_temp)) %>%
# Fit the model to each coral
mutate(ll3 = map(data, tryll3)) %>%
# Get model parameters and fitted values/residuals
mutate(pars = map(ll3, tidy),
pred = map2(ll3, data, ~augment(.x, .y)))
## Error in optim(startVec, opfct, hessian = TRUE, method = "L-BFGS-B", lower = lowerLimits, :
## non-finite finite-difference value [3]
## Error in optim(startVec, opfct, hessian = TRUE, method = "L-BFGS-B", lower = lowerLimits, :
## non-finite finite-difference value [3]
## Error in optim(startVec, opfct, hessian = TRUE, method = "L-BFGS-B", lower = lowerLimits, :
## non-finite finite-difference value [3]
pars <- initmods %>%
select(nursery, geno, pars) %>%
unnest(pars) %>%
filter(term == "ED50")
vals <- initmods %>%
select(nursery, geno, pred) %>%
unnest(pred) %>%
full_join(pars) %>%
rename(ed50 = estimate)
# Define function to plot raw data, fitted values, and ed50 for each genotype
plotfits <- function(data) {
ggplot(data = data, aes(x = max_temp)) +
geom_point(aes(y = fvfm, color = .resid > 0.15)) +
geom_line(aes(y = .fitted)) +
geom_vline(data = distinct(data, geno, ed50),
aes(xintercept = ed50),
lwd = 0.2, lty = 2) +
geom_text(data = distinct(data, geno, ed50),
aes(x = ed50, y = 0.05, label = round(ed50, 2)),
hjust = 1, nudge_x = -0.2, size = 3) +
facet_wrap(~ geno)
}
# Plot fits for each nursery
plotfits(data = filter(vals, nursery == "rrt"))
plotfits(data = filter(vals, nursery == "fwc"))
plotfits(data = filter(vals, nursery == "mote"))
plotfits(data = filter(vals, nursery == "crf"))
plotfits(data = filter(vals, nursery == "um"))
# Remove outlier datapoints
dfff <- vals %>%
filter(.resid < 0.15) %>%
select(nursery, geno, f, fm, fvfm, tank_id, tank_set, max_temp)
# Refit models individually without outliers
# Fit model to each coral, get parameters, fitted values, and residuals
fmods <- dfff %>%
nest(data = c(f, fm, fvfm, tank_id, tank_set, max_temp)) %>%
# Fit the model to each coral
mutate(ll3 = map(data, tryll3)) %>%
# Get model parameters and fitted values/residuals
mutate(pars = map(ll3, tidy),
pred = map2(ll3, data, ~augment(.x, .y)))
fpars <- fmods %>%
select(nursery, geno, pars) %>%
unnest(pars) %>%
filter(term == "ED50")
fvals <- fmods %>%
select(nursery, geno, pred) %>%
unnest(pred) %>%
full_join(fpars) %>%
rename(ed50 = estimate)
plotfits(data = filter(fvals, nursery == "rrt"))
plotfits(data = filter(fvals, nursery == "fwc"))
plotfits(data = filter(fvals, nursery == "mote"))
plotfits(data = filter(fvals, nursery == "crf"))
plotfits(data = filter(fvals, nursery == "um"))
fvals %>%
mutate(genoN = paste0("g", group_indices(dfff, nursery, geno))) %>%
ggplot(aes(x = fct_reorder(genoN, ed50), y = ed50)) +
geom_point() +
geom_errorbar(aes(ymin = ed50 - std.error, ymax = ed50 + std.error)) +
facet_grid(~ nursery, scales = "free_x") +
scale_y_continuous(limits = c(34.5, 38))
# Refit genos all together without outliers, and fixed maximum across genos
dfff$genoN <- factor(paste0("g", group_indices(dfff, nursery, geno)))
# m3f <- drm(fvfm ~ max_temp, genoN, data = dfff,
# fct = LL.3(names = c('hill', 'max', 'ED50')),
# upperl = c(100, NA, NA), # upper limit of 100 on hill parameter
# pmodels = data.frame(genoN, 1, genoN)) # all genos have fixed maximum
#
# save(m3f, file = "data/ll3modf.RData")
load("data/ll3modf.RData")
# Tidy model output and join ED50 estimates with sample metadata
res3f <- tidy(m3f) %>%
filter(term == "ED50") %>%
mutate(genoN = curve) %>%
left_join(distinct(dfff, nursery, geno, genoN)) %>%
mutate(genoN = factor(genoN, levels = genoN[order(estimate)]))
# Plot ED50 estimates for all genotypes
res3f %>%
ggplot(aes(x = genoN, y = estimate)) +
geom_point() +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error)) +
facet_grid(~ nursery, scales = "free_x") +
scale_x_discrete(labels = res3f$geno[order(res3f$genoN)]) +
scale_y_continuous(limits = c(34.5, 38))
# Get fitted/predicted values to look at fits for each genotype
pred3f <- augment(m3f, data.frame(m3f$data)) %>%
mutate(genoN = genoN.1) %>%
left_join(select(res3f, geno, genoN, nursery, ed50 = estimate)) %>% # add ed50 value for each genotype
mutate(genoN = factor(genoN, levels = levels(res3f$genoN)))
# Plot for each nursery
plotfits(data = filter(pred3f, nursery == "rrt"))
plotfits(data = filter(pred3f, nursery == "fwc"))
plotfits(data = filter(pred3f, nursery == "mote"))
plotfits(data = filter(pred3f, nursery == "crf"))
plotfits(data = filter(pred3f, nursery == "um"))
